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Abstract 

In this work we numerically study the diffusive limit of run & tumble kinetic models for cell motion 
due to chemotaxis by means of asymptotic preserving schemes. It is well-known that the diffusive limit 
of these models leads to the classical Patlak-Keller-Segel macroscopic model for chemotaxis. We will 
show that the proposed scheme is able to accurately approximate the solutions before blow-up time for 
small parameter. Moreover, the numerical results indicate that the global solutions of the kinetic models 
stabilize for long times to steady states for all the analyzed parameter range. We also generalize these 
asymptotic preserving schemes to two dimensional kinetic models in the radial case. The blow-up of 
solutions is numerically investigated in all these cases. 

1 Introduction 

Chemotaxis is one of the basic mechanisms of cell motility due to chemical interaction. Cells are attracted 
by the concentration of certain chemical substance, called chemoattractant, and they direct their movement 
toward the regions of highest concentration. Typically, this phenomenon have been described based on 
macroscopic systems of equations describing the evolution of the cell density and the chemoattractant in 
time. These systems of drift-diffusion type are the well-know classical (Patlak)-Keller-Segel models |29ll24l 
l25l . Another point of view was introduced by a mesoscopic description of these phenomena bridging from 
stochastic interacting particle systems to macroscopic equations. This middle ground consists in describing 
the movement of cells by a "run & tumble" process ll26l l28ll . The cells move along a straight line in the 
running phase and make reorientation as a reaction to the surrounding chemicals during the tumbling phase. 
This is the typical behavior that has been observed in experiments. The resulting nondimensionalized kinetic 
equation, with parabolic rescaling, reads 

?f , rr * 1 
dt 



f^ + V'VJ = - / {T e f - T*f) dv' (1.1) 



Here / = f e (t,x,v) is the density of cells at position x G M. N , moving with velocity v G V C Mr. S = S e (t,x) 
is the density of chemoattractant, whose governing equation will be described later. T e = T e [S](t,x,v,v') is 
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the turning kernel operator. We use the abbreviation /' = f(t,x,v'), T*[S](t,x,v,v') = T e [S}(t,x,v' ,v). e is 
the ratio of the mean running length between jumps to the typical observation length scale. For simplicity 
we omit the dependence on £ in notations (except T e , where £ explicitly enters into the expression.) We refer 
to ID] for the details on the rescaling. The chemoattractant S is typically given by the Poisson equation 

— AS = p where p = p e (t,x) = / f e (t,x,v)dv is the macroscopic density of cells. (1.2) 

Jv 

Here, we assume that the time scale of relaxation for the chemoattractant is much smaller that the one for the 
cell density. Although this is the most standard way of introducing the coupling with the chemoattractant, 
some authors have proposed to use some generalized models for any dimension where the chemoattractant 
S is given by 

S = - — logM*p. (1.3) 
NK 



Note that in two dimensional case, (11.31 ) is exactly (1 1 -2b - 

Starting with the kinetic equation (11.11 ). one can (at lease formally) derive the macroscopic limit as 
£ ^0, 

d tP = V_ x -(DV xP -xpV x S). (1.4) 

The details will be described in Section [3] Coupling with (1 1 -2b . one obtains the already mentioned Patlak- 
Keller-Segel system |29l l24l |25 1 for chemotaxis. We refer to the reviews lfT7ll3"01 and the references therein 
for mathematical results on this system. The behavior of the solution is quite different depending on the 
dimension. In the ID case, the global solution exists for all initial condition. In the 2D case, there exists 
a critical mass M c determined by the coefficients Ifl3l l6l. The global solution exists if the initial total mass 
satisfies M <M C (subcritical case) and its long time asymptotics are given by self-similar profiles, otherwise 
the solution blows up in infinite (critical case M = M c IH) or finite time (supercritical case M > M c lfl4l ). 
The critical case has infinitely many stationary states, each of them having their own basin of attraction H 
depending on the tail of the distribution at infinity. In the 3D case, the relevant quantity ensuring global 
existence is the I? I 2 norm of the initial density, see lfT2l . together with some condition involving the second 
moment of the solution. 

Similarly, starting with dl.lMl.3l ). one derives the modified Keller-Segel model (1 1 .4b - (l 1 -3b studied in 
[8 ], as £ — > 0. For the modified system, the existence of critical mass is extended to the ID and 3D cases. It 
is given by 

2N 2 nD 

M c = . 

X 

A very interesting consequence is that even in one dimension, this system blows-up over this critical mass. 
This makes possible to numerically analyse this highly subtle phenomena. This modified Keller-Segel sys- 
tem was numerically solved in (3l by using optimal transportation ideas and the scheme proven to be con- 
vergent in the subcritical case. 

Alt in HI 13 derived (1 1 -4b from a transport equation similar to (1 1 . lb via a stochastic model. Later the 
kinetic system dl.lb - dl.2b was formulated in [26]. Othmer and Hillen in [16, 27] studied the formal diffusion 
limit of this system by moment expansions. In [11] Chalub et al. proved that, in three dimensions with 
suitable assumptions on the turning kernel T e , the solution to the kinetic system dl.ll) - dl.2l) globally exists 
for any initial total mass. Moreover, they gave a rigorous proof of the convergence to the macroscopic limit 
for all time invervals of existence of the limiting Keller-Segel system. In the following work |[T8l |T9l the 
results are extended to more general cases. 
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The first task of this work is to design an Asymptotic-Preserving (AP) method for the kinetic system 
dl.lMl.3l > in the one dimensional case. The "Asymptotic Preserving" scheme, introduced by [20], is a 
suitable numerical method for the kinetic equation in such a way that letting e — > with the mesh size and 
time step fixed, the scheme becomes a scheme for the limiting system. In this case, this means as £ — > 0, the 
moments of the solution to the kinetic system dl.lt - jl.3b solve the modified Keller-Segel system dl.4M1.3t 
automatically. We refer to ll2~Tl and the references therein for a detailed review on AP scheme. 

As mentioned above, the solution to the modified Keller-Segel system dl.4M1.3t over the critical mass 
blows up in finite time t*. But the global solution to the kinetic system dl.lb - dl.3b exists. An attractive 
idea is that, the solution of the kinetic system after t* might be a good candidate to describe the behavior of 
solution to the Keller-Segel system after blow-up. There are other mechanisms of regularizing the Patlak- 
Keller-Segel system after blow-up by saturating the effect of the chemoattractant, introducing volume effects 
or cross-diffusion in the macroscopic system, see IjlOl for a discussion. However, to regularize it by the 
kinetic system is more appealing since the kinetic system is directly derived from interacting particle systems 
models. 

Bournaveas and Calvez [7 ] studied the kinetic system dl.lt - dl.2t with the local turning kernel T e defined 
below in (12.6b . They showed that, in the two dimensional spherically symmetric coordinate, the solution 
exists below a critical mass m c and blows up over a critical mass M c . Two questions remains unsolved. First, 
the large-time behavior of the solutions for subcritical cases is not clear yet. Second, the two thresholds 
do not match, i.e. m c < M c . Moreover, the appearance of blow-up in the one dimensional case has not 
been clarified yet, see [31]. The second aim of this work is to investigate these questions through numerical 
simulation. 

The outline of the paper is as follows. In section [2] the kinetic models we will work on are described. 
Section [3] gives the macroscopic limit of each model. Section [4] describes the AP schemes for these models. 
Finally, in section [5] we present the simulation experiments with our schemes and draw some conclusions. 

2 The kinetic models 

In this section we briefly describe the kinetic models we used in our simulations. The chemoattractant S is 
always given by the log kernel convolution (1 1 -3b . 

The turning kernel T e in the kinetic equation dl.lb needs to be specified. The turning kernel, T e = 
T e [S](t,x,v,v') > 0, measures the probability of velocity jump of cells from v' to v. To derive the Keller- 
Segel equation (1 1 .41) as e — > 0, one has to incorporate both 0(1) and 0(e) scale into T e . In the following 
work, we consider T e in the form 



Here F(y) > gives the equilibrium of velocity distribution when there is no directional preference. It is 
natural to assume 



T\ characterizes the directional preference. We can assume T\ > since we are considering the reaction of 
cells to chemoattractant. The cells have a larger probability to jump to a preferred direction. 



T e [S](t,x,v,v 



)=F(v) + eT l + 0(e 2 ). 




(2.1) 
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2.1 ID nonlocal Model 

Now we employ the nonlinear kernel introduced in ifTTI . 

T e = T e [S](t,x,v,v') = a + \if{S{t,x),S{t,x + ev)) + a-^r{S(t } x),S{t,x-ev')). 

The first term means the cell decides a new direction to move based on the detection of current environ- 
ment and probable new location, while the second term gives the influence of past memory on the choice of 
the new direction. For simplicity we neglect the second term. We take a + = 1 , a = 0, and we consider 

W(S,S)=F(v) + (S-S) + . (2.2) 

Again F = F(\v\) is the equilibrium function satisfying (12.ll ). For simplicity we introduce the notation 

8 e S(x,v) = (S(x + ev)-S(x)) + . 

Note that 8 e S = 0(e). Then the kinetic equation (11.11) reads 

£ §7 + V |£ = I { {F ( V ) + 5£ ^))P- (i + / S e S(*,v')dv')/) . (2.3) 
where i£(l = [— x max ,.x max ], v G V = [— v max , v max ], |V| = Vol(V). We impose the initial conditions, 

f(0,x,v,)=f(x,v)>0, (2.4a) 
S(0,x) =S' (x) >0, (2.4b) 

and reflection boundary condition for f, Neumann boundary condition for 5, 

/(?,±x max ,v) = f(t,±x max ,-v) , (2.5a) 
^|, =± , max =0. (2.5b) 

Equations (I2.3I )- (I1.3I ) gives the nonlocal model in one dimension. The global solution exists, regardless of 
the total mass. 

2.2 ID local Model 

In this section we summarize a turning kernel introduced by Bournaveas and Calvez Q. Let 

T e = T e [S](t,x,v,v')=F( V ) + e{vVS(x)) + , (2.6) 

where F = F(\v\) is the equilibrium function satisfying (12- 1 b - Then the kinetic equation (11.11 ) in one dimen- 
sion reads 

£ §F + v Jx~ = \ ((F(v) + £ (v ' V5) +^ 9 " (1 +Cl£ |V5|)/ ^ ' (2 - 7) 

with c\ = / (v • VS /\VS\) + dv = - / |v| dv. The same initial conditions (I2.4ab (I2.4bl ). and boundary con- 

, -/v. . . 2 7y 

ditions (I2.5ab . (I2.5bl ) are applied. 

Remark 2.1. T/'we expand the kernel ( 12.21 ) T = To + £T\ + (£ 2 ), cmd drop the terms higher than first order, 
we get exactly the kernel ( 12. 61 ). 



2.3 2D local Model 

Plug the local kernel ( 12.61 ) into the kinetic equation dl.ll ) in 2D, one obtains, 

f ed t f+vV x f=-(pF(v)-fs + ep(vVS) + -ec 2 \VS\f), v e V = B(0,v max ), x e M 2 (2g) 
1 -V 2 5 = p 

with C2 = Jy (v • VS/ 1 V5|) + dv. The equilibrium F(v) satisfies (I2.ll ). 

Now we assume the initial data f is spherically symmetric: f ! (&x, &v) = f(x,v), for any rotation 0. 
Then the solution / to (12.81 ) remains spherically symmetric. Denote 

r= \x\ e [0,°°), CO = Ivl G [0, v max ), = cos -1 — E [0,7rl. 

rco 

Then / is a function of r, ft), and Let 

h(t,r,co,d) = r/(f,x|( rj0 ),v|( fflcose£Bsine )) = r/(f,x|( rcos0 , rs ine) > v l (<» ,0)): 

p(t,r)=rp(t,r). 

Then the density is given by 

p(t,r) = rp(t,r) = r f(t,\x\=r,v)dv = 2 cohdOdco, 

Jv J VO<G)<v raax ,0<e<?r 



and the total mass is 



M = 2k [ p(t,r)dr = 4k [[( (ohdddcodr. 

JO J J ■J0<m<v mdx fl<9<7Zfl<r 



Then equation (12.81 ) can be rewritten as 

dth+j (d r (cos9h)-d e \^^ h X) = ^{pF {(o) - h) + ^{(op{cos6d r S) + - c 2 \d r S\h) , ^ ^ 
k -d r {rd r S) =p, 
where 

c 2 = 2/ / w 2 (cos 05,5/ |5 r 5|) + dwde = / / ft) 2 |cos 0| dcodB = — 
jo jo 1 Jo Jo 3 

and F(co) satisfies, 

/^'max 

271 / (oF((o)dco = 1. 

J0 

From the second equation of ( 12.91 ), d r S can be computed by 

1 r 

5 r 5 = — / pdr. 
r Jo 

We impose the boundary conditions, 

h(0,(O,d) =h(0,(O,n-d), 
dMr maxi (0,d) = 0. 

It has been shown that the solution blows up in finite time in [7j for M large enough while existing globally 
for small enough mass. 
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3 The macroscopic limits 

3.1 ID model: local and nonlocal 

The local kinetic model ( 12.71 ) and the nonlocal kinetic model ( 12.31 ) give the same asymptotic limit as e — > 0. 
We apply the Hilbert expansion into (12.71) and (12.3b . and match terms of the same order in e. The classical 

Keller-Segel system can be derived for p(t,x) = / f(t,x,v)dv, 

Jv 

d ; p = d x (Dd xP -xpd x S) (3.1) 

where 

D= [ \v\ 2 F(v)dv, and £ = i/| v | 2 dv. (3.2) 
Jv 2 Jv 

We refer to [11] for the details. Besides, by taking the moments of the kinetic equation (11.11) . one has 

^ +v -G// /dv )=°- 

Compare with (13 - lb . one has 

Dd x p-%pd x S = — / vfdv. 

e Jv 

While the reflection boundary condition (I2.5ab leads to 

/ v/dv = 0, atx = ±x max . 
Jv 

One arrives at the general boundary condition (for mass preservation) for Keller-Segel model 

DdxP ~ %pd x S = 0, at x = ±x max . 
Furthermore, the Neumann boundary condition for p is derived under the condition d2.5b| ). 

3.2 2D spherical symmetric local model 

For the reduced spherical symmetric kinetic system ( 12.91 ), one can derive the similar asymptotic limit for p 
as £ — > 0, 

d r p=d r (Drd£-xpd r S^j (3.3) 

with 

—d r (rd r S) = p 

where r G [0, r max ] with 

D = 2 / (o 3 cos 2 d F (co) dco dd = n co 3 F(co)d(o 
Jo Jo Jo 

and 

X= / w 3 cos 2 0dwd0 = ^i. 

jo jo 8 

The Neumann boundary condition is derived, 

d r S = 0, atr = 0,r max , 
- p 

d,— =0, atr = r max . 
r 
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4 The Numerical Scheme 



In this section we present the numerical method for the kinetic system (12.3l )- (12.5bl ) by the even and odd 
parity formalism, which has been successfully applied to the diffusive limit of linear transport equations in 



4.1 Odd and Even Parity 

We introduce the operator 

R[f}(x,v)=^(f(x,v)+f(x,-v)), 

J[f](x,v) = ±(f(x,v)-f(x-v)). 
The new functions R[f] and J[f] are defined in Q. x V + x R + , where V + = {v 6 V\v > 0}. Then we take 

r(*,v) = R[f] = \ (f(x,v)+f(x, -v)) , j(x,v) = J[f] = JL (/(*, v ) _ /(x , _ v )) . 



We can recover / from r and j, 

f(x,v) 



r(x,v) + ej(x,v), if v > 0, 

r(x, —v) — e j(x, — v) , if v < 0. 



4.1.1 ID nonlocal model 

Now we describe the odd and even decomposition for the nonlocal model. Since f(x, — v) satisfies 

e 9 -^ - v d Jk^ = I ((F(v) + 8'S{x, -v)) p - (l + / 8'SW) dv') /(*, -v)) , 
we obtain, 

dk+v3J = 4,((F(v)+/?[S e S])p-(l+<S e S>)r), 

d t j+-^vd x r = l(J[8 £ S]p-(l+<8 E S>)j), 

where < 8 e S >= / 5 e S(x,v')dv', and 

p=ffdv=f (f(x,v)+f(x,-v)))dv = 2 f r(x,v)dv. (4.1) 
Jv JV+ JV+ 

We assume e < 1 , then we can rewrite the equations for r and _/ as, 



1 



d t r+vd x j = ^ ((F(v) +/?[5 e S]) p - (1+ < S e S >) r) . 



a f j + v^r = 1 (/[5 e S]p - (1+ < S e S >) ;) + f 1 - -i ) vd x r. 



(4.2) 



Let us finish with the boundary conditions for new variables r and j. From the reflection boundary 
condition (I2.5ab and the definition of r and j, we can easily get the boundary condition for j, 

j(t,x b ,v) = — (f(t,x b ,v)-f(t,x b ,-v)) = 0. (4.3) 
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where xy = ±x max . Then plug this into the second equation of (14.21) . we derive , 

vd x r\ x=Xb = J[8 £ S]p\ x = Xb = ^(8 £ S(x b ,v) - 8 e S(x h ,-v))p . 

Now from the Neumann boundary condition (12.5bl ). we have S(xb + £v) = S(xb — £v), i.e. 8 e S(xb,v) = 
8 £ S(xt,,—v). So the boundary condition for r is 

d x r\ x=Xb = 0. (4.4) 

Now, we can propose an asymptotic preserving method for the one dimensional nonlocal model. The 
idea can be applied to the other models straightforwardly. 

As £ approaches 0, we can derive the Keller-Segel equation (11.41) asymptotically from the kinetic equa- 
tion (12.3b - So a natural requirement is the numerical method for (12.31 ) should discretize its macroscopic limit 
as e — > 0. We give an AP method following l22l . 

We can employ an operator splitting method on (14.2b - First the (stiff) source part is solved by the implicit 
Euler method, 

d t r = 1 ((F(v) +R[S e S])p - (1+ < 8 £ S >) r) , 

1 / 1\ (4 - 5) 

dj = -L (J[8 £ S}p - (1+ < 8 £ S >)j) + (l - i) v3,r. 

Then the transport part can be solved by an explicit method (e.g. upwind scheme), 

d t r + vd x i = 0, 
d t j + vd x r = 0. 



We can check that the method described above is AP easily. As £ — > the leading term in £ of ( 14.51 ) gives, 

F(v)+R[8 £ S] 



1+<8 £ S> 



-p=pF(v) + 0(e) 



J\8 £ S]p-vd x r / 1 . , 



Plug into the first equation of (14.61 ) and integrate over V + , we get exactly the Keller-Segel equation (11.4 
with D and # given by (13.2b . 



4.1.2 ID local model 

For the one dimensional local model, we obtain, 

d t r + vd x j = ^((F(v) + |KS|)p-(l+ Cl £|d v S|)r 

d t j + -^vd x r = ^ Qvd*Sp ~ (l+ci£ . 

The remaining work is similar to the nonlocal model. The boundary conditions for for new variables r and 
j are also given by (14.41 ) and (14.31) . The AP scheme can be designed similarly. 
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4.1.3 2D spherical symmetric local model 

Denote 

1 1 

R = -(h(r,co,d)+h(r,co,7i -d)), J = —(h(r,CQ,0) - h(r,C0,7V - 6)). 

R and J are functions defined on [0, r max ] x [0, v max ] x [0, 7l] . Then we can write the equations for R and /, 



d t R + (0 ( d,.(cos0/) -d e ( ^-^/j j = -^(pF -R) + ^ f^cop cos G\d r S\- c 2 \d r S\R 



(d r (cos6R)-dg (—^— R jj = ~^2 J+ \ (^^P co& ^ d rS-c 2 \d r S\J 



Now we split this into two steps. First the collision step can be solved implicitly, 



d t R= (pF-R) + ^ (^cop cos 6\d r S\-c 2 \d r S\R ) , 



d t J= ~^2 J +^ (^cop cos dd r S-c 2 \d r S\A + M - ) CO ( d r (cos6R) - d e (^^-R 



Then the transport part can be solved explicitly, 



d t R + CO ( d r (cos 07) - d e ( ^-J ) ) = 0, 



d t J + CO ( d r (cos 6R) - d e ( ^-R ) ) = 0. 



(4.7) 



One can show p solves the macroscopic limit (13.31 ) as e — > 0. 

4.2 Time and space discretization 
4.2.1 ID system: nonlocal and local 

In this section we give full discretization for the ID nonlocal kinetic system. The discretization for the ID 
local system is similar. First we solve (14.51) by implicit Euler method, 



^f- = ^ ((F(v)+R[8 £ S*])p* - (1+ < 8 e S* >>*) , 

= J_ {J[ 8^S*]p* - (1+ < 8 E S* >)/) + (l - 1) vdj c) r*. 

(c) 

where r n and j" are the numerical solutions of r and j at time t n . d x is the central difference discretization 
of d x . And a linear interpolation is used to get S(x + ev) and S(x — ev). < 8 e S > is evaluated by applying 
the trapezoidal rule on the interpolated value. 

If we integrate (14.51) over V + , we can get d,p = 0. So 

p* = p", S*=S". 
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Therefore r* and j* can be solved explicitly, 

» _ s 2 r n + At{F{v)+R[8 e S n })p n 
r ~ e 2 + At{\+ < 8 e S n >) ' 

£ 2 /' + At (V[S £ S"]p" + (£ 2 - l)vdj c) r*) 

7 ~ £ 2 + A?(l+<5 e S" >) 

Then we apply a first-order upwind scheme on (14.61) to get r" +1 , 



r «+i_ r « v v 

j n+l -j n v v 

+ 2Ax (r '' +1 " r, '- l) = lA^c " 2 * + h - y) ■ 

As £ — ^ 0, a simple computation shows that this scheme leads to (after integration over V + ), 
pn+ _ pn ^ . vn n^(c) Q n\ , a^/t/N^( c ) r.n 



At 



dj c) (Dd^p n -Zp n di c) S n ) + AxC(V)d£p* 



with D and % given by (13 -2b . p given by (14.11) . Here dj& is the general three-point central difference of 
d xx . Ciy) is a constant which only depends on the velocity space V. In the case of V = [— v max , v max ], 
C{V) = v max /4. 

4.2.2 2D spherical symmetric local model 

We describe a first order discretization for the transport part (14.71 ) of spherical symmetric system. By intro- 
ducing 

P{x,v) = l -{R(x,v) +J(x,v)), Q{x,v) = l -{R{x,v)-J(x,v)), 



one obtains 



d t P + <o[ d r (cosdP)-d e [ ^-P ) ) =0, 



d t Q-G)[ d r (cos6Q)-d e ( —Q ) ) =0. 



We take the grid points at 

n = (i--)Ar, i = l,...,N, 

dj = (j-^)Ad, j=l,...,M. 

The grid points and the characteristic lines are shown in FigureQ] We can define the flux at each interface 
according to the "upwind" value along the characteristic direction. An analogous idea was successfully used 
in the case of numerical simulation of the Boltzmann-Poisson system for semiconductors in [9]. Here we 
give the detailed discretization for P. The discretization for Q is similar. For each CO, 

Kt 1 -^-] ( F "+l/2J- F l 1 -l/2J G lj+l/2- G 'lj-l/2 \ _ 

At + Ar + Ad 
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e 






o r i 



Figure 1 : Illustration of grid points and characteristic lines on r — 6 plane. 



where 



with boundary 



if 7<M/2, 



(coa9j)I* +1 j, ifj>M/2+l, 
sin(0y + A0/2) sinQAe) 



'■r 



4.3 Adaptive grids for solutions with blowup 

Now let us consider the case that the initial total mass is large. It has been shown that the solutions to 
the 2D spherical symmetric kinetic model ( 12.91 ) blow up in finite time. While for the one dimensional 
kinetic local model (12.71 ). although the strict analysis is lacking iPTi . the numerical results strongly suggest 
that the solutions also blow up. In these cases, the convergence of the schemes described above would 
be questionable when simulations are performed on the fixed grids. See section for the verification of 
convergence order. 

Noting that in the blowup case, the magnitude of V X S grows up dramatically as time evolves, which 
actually characterizes another scale for the kinetic equation. To capture this scale, one has to take Ax ~ 
max l v sl • Hence the adaptive grids are needed. In numerical simulation, we double the grid points once 
||Vjf5||oo is doubled. More exactly, we apply the following algorithm. 

Algorithm 4.1. (Adaptive grids) 

s= ||V x S ||oo 

if 1 1 V X S„ 1 1 oo > Is then 

S ^ 1 1 V x^n 1 1 oo- 



n 



Double the grids, half the time step size At. 

Derive f n on the finer grids by interpolation and continue the time evolution, 
end if 

This simple idea works well when checking the blowup property and determining the blowup time. 
A nonuniform refinement might be more efficient in this problem since the mass is concentrated at some 
separated points. However this is beyond the scope of this work and is left for future study. 



4.4 A second order scheme for ID model 

One can derive a second order scheme without much difficulty for the ID (local or nonlocal) kinetic model. 
First the Strang splitting is applied. One solves the stiff source part (14.51 ) over a time step y , then solves the 
transport part (14.61) over At, then solves the stiff source part (14.51 ) over another y. 

Here we describe the second order scheme for the nonlocal model in detail. It can be applied to the local 
model without any difficulty. The stiff part (14.51 ) can be solved exactly in time. Noting that the cell density 
p and the chemical concentration S are not changed in this step. Let 

At(l+<S e S">) 



X 



then after a time step At/ 2, 



K ' 1+ < 8 e S" > y ' 
f = Xf + (1 - X) J[5 l S " ] p" + vd^r, 

with 

e 2 ) V 2 + 1+ < 8 e S" > 9 V(l+ < $ e S n >) 2 A 

Here the computation of 8 £ S n is described in Section 1431 9% is a central difference in x direction. 
Next a second-order TVD scheme is needed to solve the transport part over a time step At, 



At 2Ax " Jl+L J '- L/ 2Ax w+1 ' l - LJ 4 V Ax J 

)T X -fi , V /.* J. \ _ V (;* r,.* , .* \ V f vAt\ 



At 

where 



^ = ^ 011111110(1 (^+1 ± ^+1 ~ r ' T ^ ± ~ r *i-l^ fi-l) ■ 



Finally we can update the density p** from r** and S** is obtained. Then the stiff part is solved over another 
At/2. 



4.5 The convolution 

Finally, we apply the FFT algorithm on the computation of convolution dl.31 ) to get S. The singularity of 
log |x| at x = makes the direct numerical integral difficult, see Jl). Noting that log \x\ belongs to L 1 , we can 
avoid this problem by taking Fourier transform first. In the numerical simulation, we will always make / to 
be compactly supported in the computational region. It is not a problem about the periodicity of boundary 
condition by extending the solution to zero in a larger interval and by computing the Fourier transform there. 
In this way we avoid any kind of aliasing. 



12 



4.6 The computation of 8 E S 



To obtain a higher accuracy in computing 8 £ S in ( 12. Q , a high order interpolation is needed to compute 
S(je + ev). Here we apply the FFT based interpolation, 

where 5 is the discrete Fourier transform of S on grids Xj. 



5 Numerical results 
5.1 ID Nonlocal Model 

The first simulations are devoted to the one dimensional nonlocal model described in section IXTI 

The following simulations are set on x G Q. = [— 1, 1], v G V = [— 1, 1]. We take iV v = 64, which can 
provide good enough accuracy for numerical simulations. The constant function in V is chosen as the 
equilibrium 

F(v) = -Uy. 

In this setting, the critical mass of blow up for the limiting Keller-Segel system is 

M c = In. 

The initial conditions in the simulation are always given by, 

P \x) = Ce-^\ f e (x,v) = p I (x)F(v), (5.1) 

where C = C(M) is a constant determined by the total mass M. 

As predicted in ifTTTl . in the supercritical case M > M c , the solutions to the kinetic system converge to 
that of the Keller-Segel system only in finite time (before blow up time). After that, the asymptotic limit 
is not valid anymore. To capture the behavior of solutions to the kinetic system after that time, one has to 
resolve the small scale e. Therefore in the simulation, we need Ax = 0(e). While in the subcritical case, as 
will be shown in the following sections, the asymptotic limit seems to be valid over any time period. One 
can take Ax independent of e, as in a typical AP scheme |[2TTl . 

As for the time step length At, the simulation results suggest that, for the sake of stability, one needs 
At = eAx/ (2v max ) for long time simulation in the supercritical case. While in the subcritical case M < M c , 
as £ — > 0, the diffusive nature of the Keller-Segel system requires At = Ax 2 /2. A general choice of At would 
be 

eAx Ax 2 



At = max 



5.1.1 Convergence order of numerical scheme 



2v m;ix 2 



In this section we test the convergence order of numerical scheme described in Section 1441 We check the 
following error, 

„ ( f \ II/a*(0 -/2Ax(0Hl (Kn , 
/2Ax(0) 1 
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log 10 (Asc) 
(a) M>M C . 



i°9 10 (Ax) 
(b) M <M C . 



Figure 2: ID nonlocal model. The convergence order of scheme described in Section l4~4l for different e. 
Left: M = 4n > M c , t = 0.0025 < t b ; Right: M = % < M c , t = 0.025. 



This can be considered as an estimation of the relative error in / norm, where fy is the numerical solution 
computed from a grid of size Ax = Xmax N Xmm . The numerical scheme is said to be k-th order if < CAx^. 

The computations are performed with N x = 100,200,400, 800, 1600, 1600, At = 

Figure|2](a) shows the convergence results at t = 0.0025 for different e, with initial data given by (15.11 ) 
and total mass M = 4% > M c . The solution to the Keller-Segel system with total mass 471 blows up at 
% rj 0.0039 (see next section). The scheme shows second order convergence (in l l norm) for supercritical 
mass before blow up time. 

Figure [2] (b) shows the convergence results at t = 0.025 for different £, with initial data given by (15.11 ) 
and subcritical mass M = n < M c . The second order convergence (in l l norm) is observed for all e. 

In conclusion, our scheme has second order convergence, uniformly in e. This is a common result for 
AP scheme, see iTTSl . 

The above simulations are performed with the transport equation (14.61 ) solved by a Lax-Wendroff method. 
The use of the second order TVD method described in Section l4~4l shows a lower, but still uniform conver- 
gence in e. 



5.1.2 Global existence and finite time blow up 

Following the proof in fTTl . one can show that the solution to the kinetic system (12.3l )- (12.5bl ) is bounded on 
[0, T], for any time T. However, the Patlak- Keller-Segel system (11.41 ) can present a blow-up phenomenon in 
finite time, see (SJ. 

Now we take the initial data (15.11 ) with supercritical mass M = 4n > M c = 2%. The blowup time % 
0.0039 from the numerical simulation. 

Figure [3] shows the global boundedness of kinetic system for different e and the finite time blowup for 
the corresponding Keller-Segel model, by drawing the time evolution of the maximum value of p e and p. 



5.1.3 AP property: Convergence in £ at a finite time interval 

As mentioned before, it has been shown in ifTTl that the solution of the kinetic system can converge to that 
of the Patlak-Keller-Segel system weakly in a finite time interval [0,t*] with t* small enough. Here we 



14 




Figure 3: ID nonlocal model. Time evolution of kinetic system and Keller-Segel model for supercritical 
mass M = 4n. 



numerically check this convergence. We use the same grid size Ax = j^qq for different s in this subsection. 

In this section, we use the notations p e and f e to for the solutions to the kinetic system, while po for the 
Keller-Segel model. 

The supercritical case is studied first. We take initial data d5.ll ) with M = 4n. Figure |4] shows the 
convergence of f e — > poF as £ — > at time t = 0.002 < %, where % ~ 0.0039 is the blow up time of the 

limiting Keller-Segel model. FigureHJa) sketches the shape of p e (x) = / f(x,v)dv. As e — > 0, p e (solutions 

Jv 

to kinetic system) approach to po, the solution of the limiting Keller-Segel system. Figure |4jb) shows that 
the convergence order of f e — > p e F is almost first order (around 0.85) in I 2 norm. 

For a subcritical case, the initial data are taken to be (15.11 ) with M = n. After a relatively long time 
t = 0.01, we obtain the very similar results. See figure [5] Now we have first order convergence in e, in / 
norm. 



5.1.4 The stationary solution of kinetic system 

In this subsection we will study the solution of kinetic system with large initial mass M = 4n at a relatively 
long time t = 0.1. After a long time, the solution stabilizes toward a stationary state. This has not been 
proved nor intuitively discussed in the literature. Figure [6] shows the function ep e {ex) at time t = 0.1. 
Figure [7] shows the function ^ e , x ' v ) as a function of v at different x. These two figures suggest the stationary 
state satisfies 



£(*,v) = 




(5.3) 



for some functions p«,(x) and /L(x,v). Therefore, let us consider the ansatz 

f(t,X,v) = £f e (t,£X,v) 

p(t,x) = ep £ (t,sx) 

S(t,x) = - — log |x| *p = 5 e (ex) + C e 
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(a) Pe(x). (b) Convergence in e. 

Figure 4: ID nonlocal model. The solutions of kinetic system and Keller-Segel system at time t = 0.002, 
before blow up time. The super-critical mass M = 471 is used. The left figure shows p e (for different e) and 
p. The right figure gives the convergence order of | \f e {x,v) — p e (x)F(v)\ (2- 




* iog 10 £ 
(a) Pe(jc). (b) Convergence in e. 

Figure 5: ID nonlocal model. The solutions of kinetic system and Keller-Segel system at time t = 0.01, 
before blow up time. The super-critical mass M = 471 is used. The left figure shows p e (for different e) and 
p. The right figure gives the convergence order of | \f e (x,v) — p e (x)F(v)\\2. 
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Figure 6: ID nonlocal model. The function ep £ (ex) for different e. For supercritical mass M = 4n at time 
t = 0.1 



where C e = ^ log £, with M the total mass. The rescaled variables satisfy the following equations 

£l % + V Tx = ( f ( v ) + ^(*,v))p- (l + J 8§(x,V)dAf 
S = — — log \x\ *p 

71 

where f> 0, i £ Q £ = [-^ , v G V = [-v max , v max ]. As t ->• oo, |£ ->. 0, we have the stationary 
solution /J should solve 

= (F(v) + <5lo(x,v)),cL- ^1 + J 8L{x,v')dv'^f (5.4) 

where is obtained from p^. Clearly F(jc,v) satisfies 

/ F(x,v)dv = l and / vF(.x,v)dv = 0. 
Jv Jv 



The second identity can be derived by integrating (15.41 ) over velocity space. These two identities are also 
verified numerically. p(x) = ep(ex) is shown in Figure[6l Some snapshots of F(x,v) are shown in Figure|7] 



5.1.5 The interaction between several peaks 

It has been shown the interaction between several peaks for the modified Keller-Segel system (11.41) in Q by 
means of optimal transportation methods. Here we will check this interaction for the kinetic system. 

Case I: Two symmetric peaks, without enough mass in each peak.- The initial condition is taken as 
the sum of two gaussian-like peaks, 

yfov) = C Qe- 8 °(- - 3 ) 2 + F(v) {5 5) 
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Fi gure 7: ID nonlocal model. The snapshots of the stationary state F(x,v) — p(e'x) ^ certain locations x. 
For supercritical mass M = An at time f = 0.1 > tj,. 
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Figure 9: ID nonlocal model. The attraction between peaks, e = 0.05. The total mass is 571. The critical 
mass is 2n. Nx = 400. 



with a suitable constant C such that the total mass is 371. We take e = 0. 1. Figure [8]shows the time evolution 
of p e (left) and ||p e ||oo (right). The solution starts with diffusion separately, and then the peaks merge. At 
beginning there is no enough mass in either peak to concentrate (^ < M c = 2k). But after they merge, the 
total mass in the new peak is large enough to form an aggregation. 

Case II: Two symmetric peaks, with enough mass in each peak.- We take the same initial setting 
( 15.51 ), with a different constant C such that the total mass is 57C. Now the mass in each peak is large enough 
to concentrate (^ > M c = 2%). In Figure |9j we can see two distinguished different phases during the time 
evolution. The first phase consists in the appearance of a metastable state where the concentration in each 
peak is formed but slowly moving toward each other. Then the two peaks merge, which increases the total 
mass into the final larger peak. It continues to aggregate the mass around it and finally reaches the stationary 
state. 

Case III: Two asymmetric peaks, with enough mass in each peak.- We take a nonsymmetric initial 
setting 

f e (x,v) = C {2.2e- m{x -°^ +2.8e- 80 ^ +0 - 3 ) 2 ) F(v) 

with a suitable constant C such that the total mass is 57T. The mass in each peak is still large enough to 
concentrate and we see that the mass tends to concentrate around the center of mass located this time much 
closer to the first peak due to asymmetry in Figure [TOl 
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time time 

Figure 10: ID nonlocal model. The attraction between peaks. £ = 0.05. The total mass is 571 with 2. 271 near 
x = 0.3, 2.87T near x = —0.3. The critical mass is 27T. Nx = 400. 



iipii. 




time 

Figure 11: ID nonlocal model. The interaction between rive peaks, e = 0.05. The total mass is 117T. The 
critical mass is 27T. Nx = 400. 



Case IV: Five unsymmetric peaks, with total mass M > 5M C .- We take a nonsymmetric initial setting, 
with five different peaks, 

fi(x, V )=ciwje-^-^ 

with w = [0.5,1.2,0.8,0.6,1], c= [0.4,-0.2,-0.6,0,0.2]. Here we pick up a suitable constant C such that 
the total mass is 117C. Now, we observe in Figure QT] the complicated dynamics of merging between the 
different peaks before forming the final peak at the center of mass and converging toward the stationary 
state. 

5.2 The ID local model: blow up in finite time 

We numerically check the open problem of the blow up property of the solution to the one dimensional 
kinetic local model (12.7b - As mentioned before, a theoretical prediction on this blow up is still lacking, see 
lOTTl . We consider the initial data given by (I5.ll ). The critical mass for corresponding Keller-Segel model is 
again 

M c = 2k. 

For the super-critical case, we take total mass M = 5n. 
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Figure 12: ID local model. The convergence order of the solutions to (12.71) at different time, e = 0.4. The 
total mass is 57T > M c . 

5.2.1 The problem in the convergence with fixed grids 

In this test e = 0.4. Different grid sizes in x are used. The number of grid points in each of the simulations 
are N x = 250,500, 1000,2000,4000,8000 respectively and the time step sizes are taken as At = eAx/v max . 
We still consider the convergence order defined as in ( I5.2I ). 

Figure [12] shows the convergence order in L 1 norm at different times Tyt = kt max / 20, where 1 < k < 20. 
The solutions show a first order convergence when t < 0.12, then the convergence order decreases as time 
evolves. After t > 0. 16 a negative convergence order is seen, which means the scheme is not convergent after 
that time. One cannot improve these convergence orders by using a finer grid. Figure [12] strongly suggests 
that the blowup happens in the solution during 0.12 < t < 0.16. Next we investigate this time period by 
using the adaptive grids proposed in section 143] 

5.2.2 The convergence with adaptive grids 

Now we use the adaptive grids proposed in section 14.31 We start the simulation with different grid sizes 
N x = 500, 1000,2000,4000. Figure [T3l shows the time evolution of ||p||oo. A nice convergence toward some 
density blow-up is observed. 

5.2.3 The convergence as e — > 

Next we study the convergence as e — > 0. We apply the adaptive grids described above and take N x = 1000 
at beginning. We take total mass to be M = 5k > M c . Since the solutions to kinetic equations also blow up in 
finite time, they converge to the solution of Keller-Segel system in a totally different way. Figure [T4l shows 
the time evolution of ||p||o° of kinetic equations with different e, as well as that Keller-Segel system. The 
way of asymptotic convergence is totally different with Figure [3] where the solutions of kinetic equations 
are globally bounded. 
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Figure 13: ID local model. The time evolution of | \p\ |oo by using adaptive grids, with different initial grids. 
The circles shows the time steps when the grids are doubled, e = 0.4. The total mass is 5n > M c . 
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(a) p(r,y). 



(b) \\P(v,y)-IL 



Figure 15: ID local model. The time evolution of p for the subcritical mass M = n < M c . N x = 2000. 
£ = 0.2. 



5.2.4 Subcritical case: long time behavior 

Now we check the long time behavior of subcritical case. Let us remark that in the case of the limiting 
Keller-Segel model, it is known that the long time asymptotics should be given by a self-similar solution 
whose profile is dictated by a stationary scaled problem, see O for instance. We now check this point for 
the ID local kinetic model by taking as total mass M = %. The simulation is performed on x 6 [—10, 10]. 
A fixed grid with N x = 2000 is used, with e = 0.2. We consider the same change of variables as for the 
Keller-Segel model, 

y = 17^' Z = \ogR(t), 



R(t) = vT+27, 



R( t y 



p(t,x) 



1 



R(t) 



(5.6) 



Here p is the density we computed from numerical simulation. 

Figure [131 a) shows the time (t) evolution of p until time z = 1.5 (t = 10) as a function of y. It strongly 
suggests that p is approaching some stationary state. We denote it by p^. In Figure [TBT b) we show the time 
evolution of the relative error ||p(r,y) — p«,||i, with p^ approximated through (15.61 ) at time z = 10. Clearly 
the solutions are approaching the stationary state, in other words, 



p(t,x) 



1 



R( t y 



Rt) 



0, 



as t 



5.3 The 2D local model: solution behavior between theoretical thresholds 

In the final simulation we test the 2D local model. As shown in Q, the solution blows up with total mass 
larger than the critical mass 
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where |V| = 7T&>max = 71 ■ And the global solution exists with the total mass lower than the other threshold, 

0.8067T n _ 
m c = = 0.806. 

Besides, as mentioned in the introduction, the critical mass for the corresponding Keller-Segel model (13.31) 
is 

Mrs = = 16. 

X 

The gap between the two estimates M c and m c should be much smaller as observed from our numerical 
simulation. It is difficult to ascertain based on numerical simulations if the kinetic system show a clear 
dichotomy as in the Keller-Segel system. 

We take r max = 2, ft^ax = 1, with the equilibrium given by 

w s 1 1 

F(ffl) = J7T-. = -. 

27T/ ( f ntt todffl % 

The grid points in each directions are N r = 1000, N a = 32, Nq = 32. The initial data are taken as 

p\r) = Cre- l5r \ h'(r,(0,d) = p'(r)F((o), 

where C = C(M) is a constant determined by the total mass M. 

Figure [16] shows the time evolution of | \p\ \^/M, for different total masses M. We take e = 1. It suggests 
that the global solutions exists for M < 17, which is much bigger than the theoretical thresholds m c . While 
for mass M > 25, the solution blows up, even if the total mass belongs to the range of masses in which there 
are no theoretical results Q. Let us comment that | \p \ \„ has a upper bound due to the limitation of grid size. 
When we use a finer grid, this upper bound would be larger. 

Finally we compute the convergence of the solutions to the kinetic system as Ar — > 0, with different total 
masses M. As shown in Figure [T7J the numerical scheme shows a fist order convergence for M < 15 (note 
that the critical mass of the Keller-Segel model is Mks = 16). While for M > 25, the numerical solutions do 
not convergence, which suggests that the solutions blow up. 
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